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Using dynamical Monte-Carlo simulations we observe the occurrence of an unexpected 
shock wave in driven diffusive systems with two conserved species of particles. This U-shock 
is microscopically sharp, but does not satisfy the usual criteria for the stability of shocks. 
Exact analysis of the large-scale hydrodynamic equations of motion reveals the presence of 
an umbilical point which we show to be responsible for this phenomenon. We prove that such 
an umbilical point is a general feature of multi-species driven diffusive systems with reflection 
symmetry of the bulk dynamics. We argue that an U-shock will occur whenever there are 
strong repulsive interactions such that the current-density relation develops a double-well 
and the umbilical point becomes isolated. 

I. INTRODUCTION 

Nonequilibrium lattice gas models of interacting particles with noisy dynamics are paradig- 
matic models of systems far from equilibrium and find a wide range of applications in biological, 
social and physical contexts 0-5|. Driving forces due to bulk fields or boundary gradients lead to 
steady state currents that invalidate the condition of detailed balance and give rise to remarkable 
features which have no equilibrium counterparts. As examples we mention boundary driven phase 
transitions, spontaneous symmetry breaking and hysteresis in one spatial dimension. Particles 
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systems with two or more conserved species exhibit particularly rich behaviour 

The coarse grained space-time evolution of bulk-driven systems is governed by two fundamental 
types of excitations: shocks, which carry discontinuities, and rarefaction waves, which are con- 
tinuous self-similar solutions of the hydrodynamic limit equations Various properties of the 
fundamental excitations like stability, speed and morphology are determined by a scalar or vec- 
tor function which relates steady macroscopic currents to average particle densities, the so-called 
current density relation. The topology of the current-density function (or surfaces, in case of sev- 
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eral species of particles) such as the number of extrema and saddle points determines qualitative 
features of the large scale dynamics and in particular the number and character of the different 
stationary phases and phase transitions that one can observe in the underlying microscopic model 
0, Q] . In this way microscopic details of local particle interactions are largely irrelevant as long as 
they produce a certain type of a current density relation. 

In this work we identify a new large-scale excitation, reminiscent of a shock wave, but which 
should be unstable according to usual shock stability criteria. Focussing on models with two 
particle species we relate its appearance to a special property of the current-density relation, the 
presence of an isolated umbilic point. We shall call therefore this new excitation an umbilic shock, 
or a U-shock. The aim of this article is to describe microscopic and macroscopic properties of 
the U-shock, and to investigate conditions for its appearance and stability. We find that such an 
excitation is not at all exotic and can generally be observed in bidirectional models with left-right 
symmetry in the hopping rates, provided that there is a sufficiently strong repulsive interaction. 

The plan of the paper is the following: In Sec. HI] we introduce our model and describe the 
U-shock microscopically, highlighting its difference from a usual shock. In Sec. IIIII we discuss 
macroscopic current density relations with an umbilic point which makes the existence of the 
U-shock possible. In Sec lIVI we prove that bidirectional models with left-right symmetry in the 
hopping rates all necessarily have an umbilic point. We finish with conclusions and perspectives. 
The Appendices contain necessary technical details. 



II. THE MODEL AND A MICROSCOPIC U-SHOCK 



Our model describes particles with repulsive hard-core interaction which hop unidirectionally 
along two chains of L sites: One chain for right-hopping particles and another chain for left-hopping 
particles. At each instant of time the system is fully described by occupation numbers n^. € {0, 1} 
(for the right movers) and m^. G {0, 1} (for the left-movers). A right-moving particle at site k can 
hop to its neighbouring site k + 1 provided it is empty, with a rate that depends on the occupancies 
at sites k^k + 1 on the adjacent chain, see Fig.l. E.g. a particle hops with rate (3 if the adjacent 
sites are both occupied, etc. For clarity of presentation and analytic simplification we shall keep 
only one rate 7 = e*^ different from others, setting all remaining rates to 1, 

a = /3 = e = 1, 7 = e'' (1) 

Then the parameter 



Q = 7-l 



(2) 
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Figure 1: Bidirectional two-chain model. For solvability,the rates must satisfy a = /3 = £ = l,7 = e'' where 
ly is the interchain interaction constant 11|. Coupling to boundary reservoirs is indicated by boxes marked 



L (the left reservoir) and R (the right reservoir). 



which ranges from —1 to oo, measures the interaction strength between the left- and right-moving 
species. For Q = the model reduces to two independently running totally asymmetric exclusion 
processes. The reason for the given choice of rates is a simplification that it offers: The current- 
density relation can be found analytically as explained in Sec. IIIII and can therefore be analyzed in 
detail. For monitoriiig the microscopic position of shocks on the lattice we also introduce a second 
class particle (SCP) [10|. This is a phantom particle which is located at some site k (not disturbing 
real particles) and hops to the right — t- A; + 1 if it finds a configuration n^+i = 1, fn-fc+i = on its 
right, or it hops to the left k ^ k — 1. It is well-known that such dynamical rules favor positioning 
of the SCP at the middle of a local density gradient that corresponds to a shock on macroscopic 
scale. 

The bulk dynamics is complemented with boundary conditions: We consider open boundaries 
where at the left end of the chain a right mover can enter the chain and it can leave it at the right 
end. Left movers are hopping to the left with the same dynamic rules. Note that in general we do 
not require complete left-right symmetry so the entrance and exit rates for different species can be 
different. We choose a maximal feeding regime where we put a particle on the entrance site once 
it becomes empty and take it out from the exit site once it reaches it. The unidirectional hopping 
along with the open boundaries ensure that a non-zero steady state current is maintained. 

Our results do not depend qualitatively on how exactly the maximal feeding regime is realized. 
For our dynamical Monte-Carlo simulation we choose the following random sequential update 
procedure. For a chain of length L, i.e. a system of 2L sites (numbered i = 1,2,...L for right 
movers and i = L + 1, L + 2, ...2L for left movers) one Monte-Carlo step consists of 2L + 2 uniform 
drawings of a random number s in the range < s < 2L + 1. If s falls within a segment [0, L], the 



4 



configuration of right movers is updated. If s = 0, and the left boundary site i = 1 is empty, we fill 
it with a particle (free entrance). If s = L and the respective site contains a particle, we remove it 
(free exit). For intermediate < s < L, if site s contains a particle, a hopping is performed on the 
right neighbouring site with given rates ([1]), provided it was empty. The update of the left movers 
is done analogously. We start from an empty lattice and after a transient time we measure site 
occupancies n^, m^, and take averages over many Monte Carlo steps and many histories. Typically 
we choose a system size up to L = 1000 sites in each chain. The transient time for L = 1000 is 10^ 
Monte Carlo steps, and averaging over up to 10 histories is done. We perform the measurements 
for different values of the interaction parameter Q. Note that due to the hardcore exclusion, the 
average densities of the right and left-moving particles may only take values between and 1. 

The maximal feeding regime usually leads to the largest particle current since we facilitate 
maximally the entrance and exit of particles at the boundary. In the absence of interaction (Q = 0) 
such boundary conditions lead to a state with average particle densities 1/2 a state with maximal 



possible particle current 
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la ]. In the presence of interaction Q the stationary density profile 



does not undergo qualitative changes for a vast interaction range —0.75 < Q < oo, see Fig 2(a) 



However, for values of Q < —0.75 one observes something very unusual and different. The bulk 
density profile becomes inhomogeneous and consists of two plateaux with an interface in the middle. 
The profiles of the two species are left-right symmetric but in each plateau the densities pi,P2 of 



the left- and right-movers are different, see Fig 2(b) As the interaction becomes stronger, the 
difference pi — p2 grows and reaches the maximum pi — P2 = 0.5 for the extreme case Q = —1. 
Note that the asymmetry of the profile is not a result of a spontaneous symmetry breaking since 
the profiles are left-right symmetric and the stationary currents of both species remain equal. 

What is the nature of the observed state? The maximal feeding regime in a bulk-driven particle 
system usually produces steady states that are controlled by rarefaction waves which (in an infinite 
system) are self-similar solutions of the type p{x,t) = p{ ^~^° )-, see also Appendix A. However, the 
interface in the middle cannot be a rarefaction wave because it does not change with time. With an 
increase of the system size, the interface gets wider. However the widening is due to a fluctuation 



of a position of the interface, as an analysis using second class particle shows, see Fig. (2(b) ). This 
is a property typical of a shock. 

However, the interface we observe is not a usual shock, either. In order to see this, it is instructive 
to look at individual particle trajectories across the interface, see Fig. ([3|). Unlike the trajectories 
running across a shock, the particles are moving slowly (in environment of large density) across 
the left side of the lattice, and then accelerate after crossing the inhomogeneity to the right-hand 
side of the lattice. Moreover, according to usual shock stability conditions, using characteristic 
velocities (see details in Sec. IIII[) . the interface shown on Fig. ([2]) should be unstable. As we shall 
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Figure 2: Average density profiles for right movers (thick hnes) and left movers (thin lines), above the phase 
transition Q = —0.4 (Panel(a)) and below the phase transition Q = —0.9 (Panel(b)). The dotted curve on 
Panel(b) shows the U-shock viewed from a second class particle. The U-shock profile, seen from the SCP 
position, does not depend on the system size (data not shown). Straight line pieces, showing bulk densities 
predicted from Eqs. dS]),®, serve as a guide for an eye. 
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Figure 3: Space-time trajectories of the right movers across a U-shock. Every 10-th trajectory is shown. A 
system of 1000 sites was equihbrated for 2 * 10^ Monte Carlo Steps before the trajectories were recorded. 
The parameters are: Q = —0.9. 

argue below, the reason for the stability and existence of the new state is an isolated umbilic point 
in the current-density relation. 



III. UMBILIC POINT IN A CURRENT-DENSITY RELATION 

The model ([T]) that we consider has the remarkable property that the stationary distribution is 
a product measure [h]. The steady state probabilities of any configuration C in a periodic system 
are given by 

L 
k=l 

where nk,mk are particle occupation number on site k on chains 1 and 2. With ([3|) the stationary 
currents ji and j2 of the right- and left-moving species can be calculated exactly as 

ji{u,v) = u{l -u) + QQii{u,v)noo{u,v) (4) 
j2iu,v) = -v{l -v) - Qnu{u,v)^looiu,v) 
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where u and v are the average densities of the right and of the left movers. The quantities Jin and 
f^oo are are stationary probabihties to have two adjacent particles and two adjacent holes, given 
by 



[u + V - l)Q - I + J {{u + V - l)Q - ly + AQuv 
= ^ (5) 

Ooo = I — U — V — Jill. 

From the stationary currents we construct the flux Jacobian (Dj) 

( dji dji\ 
du dv I 

Its two eigenvalues cx^^iu^v) play a fundamental role as characteristic speeds of the system of 
conservation laws 

dtu^d^h{n,v)=^ (7) 

which describes the coarse-grained dynamics on macroscopic scale. Microscopically the characteris- 
tic speeds are the velocities of the localized perturbations of a stationary homogeneous background 



111 ]. As such, they determine stability of shocks and rarefaction waves in the 



with densities u, v 
system. 

A commonly made assumption about the flux functions ji, j2, called strict hyperbolicity, reads: 



lyperbolic systems 



14( 1 . As argued in the 



the characteristic speeds are different ci{u,v) ^ C2{u,v) for all u,v. Strict 
have only two types of fundamental solutions: shocks and rarefaction waves 
previous section, the [/-shock is neither a usual shock nor a rarefaction wave, so it cannot not be 
stable in a strictly hyperbolic system. 

Indeed, our system is not a strictly hyperbolic one, but it has a so-called umbilic point which is 
defined as a point in the u — v density plane where where the two characteristic velocities coincide. 
It is straightforwardly verified from the analytic expressions for the currents that for our system 
this is the case at u* = v* = 1/2, where the two characteristic speeds 2 equal and zero for 
all values of Q. 

For a full discussion of the current-density relation ^ and the associated flux Jacobian ([6]) we 
note that the points where one characteristic speed vanishes generically correspond to a family 
of rarefaction waves , see also Appendix A. Looking at the location of the points where 

at least one characteristic speed vanishes, we find two different topologies, depending on Q. For 
the interaction range —3/4 < Q < 00 the umbilic point {u*,v*) is a crossing point of the curves 
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Figure 4: Location of the curves in u — w plane where at least one characteristic velocity vanishes Ci(u, v) = 
0, for small negative Q = — 0.5, — 0.75( Panel (a), thick and thin lines respectively) and large negative 
Q = —0.8, —0.9, — 0.99( Panel (b), thin, medium and thick lines respectively). The point m = t; = 1/2 is 
an umbilical point where ci = C2 = for any value of Q. For Q < —0.75, the umbilical point becomes an 
isolated point. 

The appearance of the isolated umbilic point is a consequence of a change of topology of the 
current surfaces jk{u, v, Q) from a convex to a saddle point shape at a critical value of Qcr = —3/4. 
To understand from a microscopic perspective how this happens consider a cut of the current 
surface along the line v = 1 — u for small ^ = Q + 1 <^ 1. At half- filling u = v = 1/2 the 
system comes into a configuration where all adjacent sites are either both occupied or both empty 
and gets essentially stuck: all hoppings from this configuration are suppressed by a small hopping 
rate 7^1. This strong repulsive interaction between the particles on the two lanes explains the 
occurrence of a double maximum in the curve j{u, 1 — u,Q) for Q < Qcr as shown on FiglS) The 
appearance of the double maximum then gives rise to the isolated umbilic point. The positions 
and amplitudes of the extrema of g{u,Q) = j{u,l — u,Q) are readily found from (j4]), ([5]): for 
Q > Qcr the j{u,l — u,Q) curve has one maximum at Uq = 1/2, while for —1<Q< Qcr it has 
three extrema at positions Uq = 1/2 and 2 = ^ =1= \/3(5~^ + 4/4. The respective currents are 
j* = ^/QTT/(2VQTT+ 2) and = = 1/{8\Q\). 

Now we are in a position to analyze the U-shock solution. We can identify the average bulk 
densities of the left plateau of the U-shock (see Fig. 2) with ul,U2 respectively, for the 

following reasons: (A) We expect the U-shock plateaux, piecewise, to be governed by a rarefaction 
wave, meaning that at least one of characteristic speeds must vanish, Cj(u_, i;_) = 0. (B) We expect 
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Figure 5: Bidirectional two-chain model. Cuts of the stationary current surface along v ~ 1 — u line: 
ji{u, 1 — u,Q) at different Q = —0.5, —0.75, —0.9, —0.9999 (curves up to down). The cuts of the stationary 
current surface along the perpendicular direction u ~ v remain convex for all values of Q (data not shown) . 



the stationary currents amplitudes for right and left-moving species to be equal, due to left-right 
symmetry, i.e. ji{u-,V-) = — j2(n_,?;_). We readily find, using (jl]),([5]), three pairs of solution 
satisfying (A)and (B): (i)(n__,7;_) = (1/2,1/2); (ii) {u^,v^) = {ul,U2) and (iii) (n_,w„) = (ujjUj). 
Comparing with the U-shock, we find (m_,u_) = (nj,n2) to be the relevant solution. Indeed the 
solution (iii) is not compatible with our boundary conditions, while the solution (i) would result 
in a reduction of the particle current and is dynamically unstable. 

For the right plateau of the U-shock, we find analogously (^4., v+) = {u2, uj). Apparently, the 
first solution (i) is unstable for Q < Qcr- So, we have for Q < Qcr 

V3Q^i +4 



u^{Q) = v+{Q) = - + 



v.{Q) = u+{Q) 



The stationary currents for the U-shock are then given by the jfc(uJ,U2,Q), 



Jstat 



(8) 
(9) 

(10) 



Above the critical point, the bulk densities are u - 
by jfc(l/2,l/2,Q), 

VQ + T 



1/2, and the stationary currents are given 



Utnt = — ; for Q > . 

2 VQTT + 2 ^ - 4 



(11) 



Our analytical predictions are well borne by the MC simulations, see Figl6j Note that by 
establishing an inhomogeneous state (the U-shock) below Qcr the system optimizes its current, 
which would be strongly suppressed for any symmetric bulk homogeneous state. This can be 
viewed as a generalization of the phenomenon of a current maximization at maximal feeding [l^ 
to a system with two species. 
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Figure 6: Stationary currents of the bidirectional model at maximal feeding as function of Q. Circles show 
Monte-Carlo results for a system of size 600, while solid curve is the theoretical prediction ([T0| . ([Tl]) . Below 
the critical Q, the branch ([TIT) , indicated by the broken line, becomes unstable. pT|) . 



We can study the phase transition in our model aX Q = Qcr- Let us choose the difference 
U-{Q) — V-{Q) = A between the bulk densities of the right and left movers to be our order 



parameter. We have A = for Q > Qcr and A = ^\/3Q -"^ + 4 for Q < Qcr- Near the transition, 



^\Qcr-Sq ~ V "^9/3) so we have a square root singularity, similar to that arising in the Landau 
theory of continuous phase transitions generated by a change of a free energy potential from a 
single to a double minimum topology. On the other hand, the stationary current is continuous 
across the transition point together with its first derivative, 

jstatlq^^.Sq - jstatlq^^ =0{5q^), for (5g < 1. 

Finally, we comment on robustness of the U-shock. The U-shock turns out to be very robust 
with respect to a choice of the boundary conditions (BC). We observe the U-shock for a wide 
choice of BC, which in particular do not need to be left-right symmetric. In order to formulate the 
conditions for its appearance more precisely, we need to parametrize the BC in terms of boundary 
reservoirs. Such a parametrization involves further technical details |171 | and is out of the scope of 
the present paper. 

The U-shock is also stable w.r. t. a change of the model parameters, e.g. the bulk hopping 
parameters, as long as they remain left-right symmetric. In particular, the particle-hole symmetry 
of the hopping rates ([T]), which makes the particle current invariant w.r.t. u,v — t- 1 — u, \ — v 
interchange (see Figl5]) can be relaxed without causing qualitative changes to a U-shock. 
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IV. UMBILIC POINT IN BIDIRECTIONAL MODELS 

With the system we have studied, we were lucky enough to find the stationary currents ana- 
lytically and establish the existence of the umbilic point. How exotic is this point? Moreover, is 
it possible to predict, from the microscopic transition rates, if the system will have such a point? 
In this section we prove that an umbilic point with ci = C2 = is not at all exotic and is present 
necessarily in bidirectional models with left-right symmetry. 

Let us consider Markov processes involving two driven particle species which are biased in 
opposite directions, their bulk dynamics (but not necessarily their boundary dynamics) being 
invariant under the left-right interchange. Let us denote the average particle densities of the two 
species as u, v and the respective currents as ji{u,v), j2{u,v). From the left-right symmetry we 
have j2{u, v) = —ji{v, u). Let us consider the line of equal densities u = v. Along this line the flux 
Jacobian (IHl) takes the form 



Dj 



a{u) b{u) 
— h{u) —a[u) 

with the respective eigenvalues (characteristic speeds) 



ci{u = v),C2{u = v) = =F\/ — 6^. (12) 

Let us assume, in addition, that we have a restriction on a number of particles which can occupy 
a single lattice site, i.e. that the maximally allowed density of the each species is limited to the 
same value maxu = maxv = Pmax- Li many applications such a restriction is a consequence of the 
hard core exclusion interaction. 

Now, let us move along the line u = v from u = to the maximally allowed value Umax and 
assume without loss of generality that our process evolves on two parallel chains, Umax = 'Wmax = 
Pmax- Let us call the lane with right movers the lane A, and the lane with left movers the lane B. 

Guided by the physical meaning of the characteristic speeds as the velocities of localized pertur- 



bations of a homogeneous state [111] , we deduce that a stationary state with small density of right 
movers on one lane and small density of left movers on another lane, attainable for u = v ^ 0, 
has characteristic speeds of opposite signs in accordance with (|12p. In this limit the left and right 
moving species are practically uncoupled and we can attribute the positive characteristic speed to 
right-moving particles ca{u = v ^ 0) > and the negative characteristic speed to the left-moving 
particles cb{u = v ^ 0) < 0. At the other end of the line u = u — )• Pmax we have vanishing density 
of left moving holes in the dense background of the right moving particles on lane A, and similarly 
for the right moving holes on the lane B. Repeating the arguments for the characteristic speeds, 
we have ca{u — )■ Umax) < and cb(u — )• Umax) = —ca{u UrasLx)- By continuity we deduce that 
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there exists a point at which the characteristic speed ca changes sign ca{u* = v*) = 0. Moreover, 
due to (jl2p . the other characteristic speed at this point, also vanishes cb{u* = v*) = 0. The point 
u* = V* where both characteristic speeds vanish is an umbilic point. 

This argument estabhshes the existence of at least one umbilic point with ci = C2 = for left- 
right symmetric bulk dynamics. Notice that the boundary conditions do not enter the argument. 
Hence one may impose boundary conditions that violate the left-right symmetry with destroying 
the umbilic point. 



V. CONCLUSIONS 



We have considered an open bidirectional two-component driven diffusive system with left-right 
symmetry for the bulk dynamics in the maximal flow regime. We have discovered and described 
in detail a bulk inhomogeneous solution of a novel type, denoted U-shock. This solution has many 
properties of a usual shock - in particular its microscopic sharpness - but does not satisfy the 
usual criteria for the stability of shocks. We have computed the critical value of the (repulsive) 
interaction above which the U-shock exists. We have shown that the existence of the U-shock is 
due to the intrinsic presence of an isolated umbilic point with vanishing characteristic velocities. 
No fine-tuning of the interaction strength is required above the critical value. 

The U-shock turns out to be robust also with respect to changes in the boundary parameters, 
even if they violate the (necessary) left-right symmetry of the bulk dynamics. In an open system 
it is only necessary maintain a stationary a maximal flow regime. More generally, we proved the 
existence of an umbilic point with vanishing characteristic velocities in any bidirectional model 
with left-right symmetry of the bulk hopping rates. The presence of umbilic points alters crucially 
the dynamics of the system in the hydrodynamic limit ()Aip . givi ng rise to a variety of unusual 
solutions, called undercompressed and overcompressed shocks 18|. The U-shock is one of such 



solutions. The necessary condition for an U-shock is a sufficiently strong repulsive interaction, 
resulting in a saddle point in the current-density surfaces. 

Bidirectional models are being widely studied in the literature, in particular, in connection with 
the intriguing phenomenon of spontaneous symmetry breaking (SSB) jl9t|-[28l|. However, we find 
that in most known cases the current-density relations are convex surfaces. It would be interesting 
to study SSB in presence of an isolated umbilic point. Another interesting problem is to explore 
the full phase diagram of the open system and to find out which role the U-shock might play in 
boundary driven phase transitions. 
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Appendix A: Rarefaction-wave controlled stationary states 

The generical importance of the points where at least one of characteristic speeds q vanishes 
can be demonstrated by the following argument. It is well-known [3] that partial differential 
equations of the type 

dpk . djk{pi,p2,-pK) _ „ . . -I X 

k= 1,2,.K, 

where K is the number of species, admit two fundamental classes of solutions: shock waves and 
rarefaction waves. A rarefaction wave is a self-similar solution of ()Aip . depending only on the ratio 
^ = (x — xo)/t where xq is the position of its center, and t > 0. We argue that in the long-time 
(stationary) limit t — t- oo the stationary bulk density pstat generated by a rarefaction wave has zero 
characteristic speed Cp{pstat) = 0. Here p{x,t) is a vector the components of which are the density 
profiles pi{x,t), p2{x,t), pK{x,t) of the respective species. By pstat we denote the vector with 
stationary bulk densities {pf '^^ /o|*"*, pf^"*}. 

We search for a solution of (lAlj) in the form p{x,t) = /i(^). Substituting in (lAlj) . we obtain 

where the matrix (Dj) is the Jacobian of the flux (-Dj)pg = djp/dpq. The above equation can be 
rewritten as 

In the limit i — )• oo, the scaled displacement ^ = (x—xo)/t — > vanishes for any finite x — xq, and 
the above equation reduces to iDj)\^^^ h' = 0. Each solution of this equation is an eigenvector 
of the flux Jacobian Dj with zero eigenvalue. Consequently, the matrix (-Di)t^oo = {D]){pstat) 
is a matrix with zero eigenvalue, i.e. at the point pf,tat at least one Cp{pstat) = 0. The respective 



rarefaction wave is called p-rarefaction wave 



151 ]. Of course, in order to guarantee the stability 



of the above discussed rarefaction wave with respect to local perturbations at the boundaries, the 
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boundary conditions have to be chosen appropriately. 
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